Probing the relaxation towards equilibrium in an isolated strongly correlated ID Bose gas 
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The problem of how complex quantum systems eventually come to rest lies at the heart of statistical me- 
chanics. The maximum entropy principle put forward in 1957 by E. T. Jaynes suggests what quantum states 
one should expect in equilibrium but does not hint as to how closed quantum many-body systems dynamically 
equilibrate. A number of theoretical and numerical studies accumulate evidence that under specific conditions 
quantum many-body models can relax to a situation that locally or with respect to certain observables appears 
as if the entire system had relaxed to a maximum entropy state. In this work, we report the experimental obser- 
vation of the non-equilibrium dynamics of a density wave of ultracold bosonic atoms in an optical lattice in the 
regime of strong correlations. Using an optical superlattice, we are able to prepare the system in a well-known 
initial state with high fidelity. We then follow the dynamical evolution of the system in terms of quasi-local den- 
sities, currents, and coherences. Numerical studies based on the time-dependent density-matrix renormalization 
group method are in an excellent quantitative agreement with the experimental data. For very long times, all 
three local observables show a fast relaxation to equilibrium values compatible with those expected for a global 
maximum entropy state. We find this relaxation of the quasi-local densities and currents to initially follow a 
power-law with an exponent being significantly larger than for free or hardcore bosons. For intermediate times 
the system fulfills the promise of being a dynamical quantum simulator, in that the controlled dynamics runs for 
longer times than present classical algorithms based on matrix product states can efficiently keep track of. 



Ultracold atoms in optical lattices provide highly control- 
lable quantum systems allowing to experimentally probe var- 
ious quantum many-body phenomena. In this way, ground 
state properties of Hamiltonians that play a fundamental role 
in the condensed matter context have been investigated un- 
der precisely tunable conditions |[I}j3|. Features that are even 
harder to probe in actual condensed matter materials or to 
simulate in numerical studies are dynamical ones, including 
dynamical properties emerging in adiabatic sweeps 13 and 
far from equilibrium ifS HTTIl . In this respect, for example, the 
quench from a shallow to a deep optical lattice |6-8| and the 
phase dynamics emerging after splitting a one-dimensional 
Bose liquid 1 12 1 have previously been studied experimentally. 

In this article, we report on the direct observation of re- 
laxation dynamics in an interacting many-body system using 
ultracold atoms in an optical lattice. Starting with a patterned 
density with alternating empty and occupied sites in isolated 
Hubbard chains, we suddenly switched on the tunnel coupling 
along these chains and measured the emerging dynamics in 
terms of quasi-local densities, currents and coherences. Both 
the initial state preparation and the detection was realized us- 
ing a bichromatic optical superlattice |[T3l fT4l . For a wide 
range of (repulsive) inter-particle interactions, we find a fast 
relaxation of the measured observables to steady state val- 
ues which are consistent with a dynamical version of Jaynes' 
principle 115]. The timescale of the relaxation cannot be at- 
tributed to a classical ensemble average. For short times, we 
compare the experimental results to time-dependent density- 
matrix renormalization group simulations (t-DMRG, for a re- 
view see Refs. 1 16, 17 1 and references therein) of the Hamil- 
tonian dynamics without free parameters, further developing 
the ideas of previous numerical studies U8ul9i . 



Concept of the experiments. We consider a one- 
dimensional chain of lattice sites coupled by a tunnel coupling 
J and filled with repulsively interacting bosonic particles. In 
the tight-binding approximation, the Hamiltonian takes the 
form of a one-dimensional Bose-Hubbard model [3, 20J 
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where dj annihilates a particle on site j, hj — a}jdj reflects 
the number of atoms on site j and U is the on-site interaction 
energy. The parameter K — muj'^(P (m is the particle mass, 
d the lattice spacing) describes an external harmonic trap with 
trapping frequency w ~ 27r x 61 Hz, present in the experi- 
ments. 

The experimental sequence can be described in three parts 
(see Fig. la): (/) At t = 0, the system is initialized in a 
density wave represented as a state vector \tlj{t = 0)) — 
I • • • , 1, 0, 1, 0, 1, • • • ), such that only lattice sites with an 
even site index are occupied and no tunnel-coupling is present 
along the chain. (//) After the quench to a distinct set of pos- 
itive parameters J, U and K, the system follows the non- 
equilibrium dynamics of the above Hamiltonian H. (Hi) Fi- 
nally, the tunnel-coupling is suppressed again and the proper- 
ties of the evolved state vector \ip{t)) are read out. 

We started our experiments by loading a BEC of about 
45 X lO'^ ^""Rb atoms in the \F — 1, mp — —1) Zeeman level 
into a 3D optical lattice formed by retroreflected laser beams 
of wavelength Xxi = 1530 nm along one direction ("long lat- 
tice") and Xy^z = 844 nm along the other two. In this loading 
we were crossing the transition to a Mott-insulator which re- 
sulted in an occupation of not more than one particle per site. 
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FIG. 1. Relaxation of the density pattern, (a) Concept of the exper- 
iment: after having prepared the density wave \tp{t = 0)) (;), the 
lattice depth was rapidly reduced to enable tunneling (/(). Finally, 
the properties of the evolved state were read out after all tunneling 
was suppressed again (Hi), (b) Even-odd resolved detection: parti- 
cles on sites with odd index were brought to a higher Bloch band. A 
subsequent band-mapping sequence was used to reveal the odd- and 
even-site populations 1 13, 14 1. (c) Integrated band-mapping profiles 
versus relaxation time t for h/{AJ) ~ 0.9 ms, U / J = 5.16(7) and 
K/ J ~ 9 X 10^'^. (d) Odd-site density extracted from the raw data 
shown in c. The shaded area marks the envelope for free Bosons 
(light grey) and including inhomogeneities of the Hubbard parame- 
ters in the experimental system (dark grey). 



Finally we added to the long lattice another optical lattice with 
wavelength X^g = 765 nm = Xxi/"^ ("short lattice") with the 
relative phase between the two adjusted to load every second 
site of the short lattice 1,14. 21 1. Completely removing the 
long lattice gave an array of practically isolated ID density 
waves I^at) = | • • • , 1, 0, 1, 0, 1, • • • ) - thus realizing step (/) 
- with a distribution of particle numbers N and thus lengths 
L — 2N — 1 given by the external confinement. For our pa- 
rameters, we expect chains with a maximal particle number of 
^max — 43 and a mean value of iV ~ 31 (see Supplementary 
Material for details on the loading procedure). 

To initialize the many -body relaxation dynamics of step (//), 
we quenched the short-lattice depth to a small value within 
200 /IS, allowing the atoms to tunnel along the x-direction. 
After a time t, we rapidly ramped up the short lattice to its 
original depth, thus suppressing all tunneling. Finally, we 
read out the properties of the evolved state in terms of den- 
sities, currents and coherences in step (///). Note that in the 
experiments we always measured the full ensemble average 
X{t) — E[N}{ipN{t)\X\ipN{t)) of an observable X over the 
array of chains (denoted by the averaging operator E{jv}), 
rather than the expectation value for a single chain with N 
particles. 

Relaxation of quasi-local densities. We first discuss mea- 
surements of the density on sites with either even or odd index. 
After the time evolution, we transfeiTed the population on odd 
sites to a higher Bloch band using the superlattice and detected 
these excitations employing a band-mapping technique (see 
Fig. lb) lfT3l[T4 1. Fig. Ic shows the integrated band-mapping 
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FIG. 2. Relaxation of the local density for different interaction 
strengths. We plot the measured traces of the odd-site population 
"-odd(i) for four different interaction strengths U/J (circles). The 
solid lines are ensemble-averaged results from t-DMRG simulations 
without free parameters. The dashed lines represent simulations in- 
cluding next-nearest neighbor hopping with a coupling matrix ele- 
ment Jnnn/J ==; 0.12 (a), 0.08 (b), 0.05 (c) and 0.03 (d) calculated 
from the single-particle band structure. 



profiles as a function of relaxation time for h/{AJ) ~ 0.9 ms, 
U/J = 5.16(7) and K/J ~ 9 x IQ-^. We plot the resulting 
traces rioddlO ™ Fig- 1^. We generally observe oscillations 
in riodd with a period T ^ h/ (4 J) which rapidly dampen out 
within 3-4 periods to a steady value of ~ 0.5. The same qual- 
itative behavior is found in a wide range of interactions (see 
Fig. 2). 

We performed i-DMRG calculations, keeping up to 5000 
states in the matrix-product state simulations (solid lines in 
Fig. 2). The Bose-Hubbard parameters used in these sim- 
ulations were obtained from the respective set of experi- 
mental control parameters. Furthermore, we took into ac- 
count the geometry of the experimental setup by perform- 
ing the corresponding ensemble average E{jv} over chains 
with different particle numbers N (see Supplementary Ma- 
terial). For the times accessible in the simulations, these av- 
erages differ only slightly from the traces obtained for a sin- 
gle chain with the maximal particle number A^max = 43 of 
the ensemble (see Supplementary Material). For interaction 
strengths [// J < 6 (Fig. 2a-c), we find a good agreement 
of the experimental data and the simulations. In this regime, 
only small systematic deviations can be observed, which are 
strongest for the smallest value of U/J which coiTesponds 
to the smallest lattice depth. They can be attributed to the 
breakdown of the tight-binding approximation for shallow lat- 
tices which gives rise to a significant amount of longer-ranged 
hopping. When including a next-nearest neighbor hopping 

term — Jnnn X]j('^]'^j+2 + h-c) in the t-DMRG simulations 
we obtain quantitative agreement with the experimental data 
(dashed line in Fig. 2). For larger values of U/J and coiTe- 
spondingly deeper lattices, the tight-binding approximation is 
valid. For U/ J > 10 (Fig. 2d), larger deviations are found. 
Here, the dynamics become more and more affected by resid- 
ual inter-chain tunneling and non-adiabatic heating as the ab- 



solute timescale of the intra-chain tunneling cx 1/ J becomes 
larger. 

The results of the density measurements can be related to 
the expectations for an infinite chain with K — 0. There, 
the time-evolution can be calculated analytically in the case 
of either non-interacting bosons ([// J = 0) or infinite inter- 
actions (U/J — > oo) ITSl fl9l . These limiting cases can be 
well-understood through the mechanism of local relaxation by 
ballistically propagating excitations. The on-site densities fol- 
low 0-th order Bessel functions describing oscillations which 
are asymptotically dampened by a power law with exponent 
— 1/2. The damping we observe in the interacting system, 
however, is much faster This behavior has also been found in 
t-DMRG simulations of homogeneous Hubbard chains with 
finite interactions |18 19|. The exact origin of this enhanced 
relaxation in the presence of strong correlations constitutes 
one of the major open problems posed by the results presented 
here. 

Measurements of quasi-local currents. Employing the 
bichromatic superlattice, we were also able to detect the mag- 
nitude and direction of quasi-local density currents. Instead of 
raising the short lattice at the end of step (//), we ramped up 
the long lattice to suppress the tunnel-couplings through ev- 
ery second potential barrier in the chain (see Fig. 3a). At the 
same time, we set the short lattice to a fixed value to obtain 
always the same value of {U/J)uw — 0.2 in the emerging 
double wells. By tuning the relative phase between the long 
and short lattice we were able to selectively couple sites with 
index (2j, 2j+l) ("even-odd", j integer) or (2j— 1, 2j) ("odd- 
even"). We recorded the time-evolution in the now isolated 
double- wells using the same final read-out scheme as for the 
densities (see Fig. 3b). We find sinusoidal tunnel oscillations 
which dephase only slowly and decrease in amplitude with in- 
creasing relaxation time t. The phase 4> and amplitude A of 
these oscillations were extracted from a fit of a sine-wave to 
the data and are plotted in Fig. 3c as a function of the relax- 
ation time for U/J = 5.16(7). While the phase contains the 
information about the direction of the mass flow, the ampli- 
tude is a combination of the local population imbalance and 
the strength of the local current. 

We find cf) to evolve linearly in time, giving strong evidence 
that the excitations in the system expand approximately ballis- 
tically as suggested in Refs. liT8l[T9l . Furthermore, its value 
does not change when coupling even-odd or odd-even sites, 
indicating the absence of center-of-mass motion in the sys- 
tem. The amplitude A on the other hand decays to zero on 
the same timescale as the oscillations in the local densities 
dampen out - in fact the quantities (1 ± A)/2 provide en- 
velopes to the traces riodd and novcn (see Supplementary Ma- 
terial). On short timescales, < AJt/h < 3, we find the 
decay of the amplitude to follow an approximate power-law 
cx with a — 0.86(7). This behavior might well change at 
longer times, where no significant amplitude was measurable. 
We extract the power-law coefficients a for a wide range of 
U/J (right inset to Fig. 3c). In all cases, the absolute values 
of the coefficients are larger than the one expected for free par- 
ticles, where a = 0.5, again indicating the faster relaxation in 
the presence of interactions. 
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FIG. 3. Quasi-local current measurement, (a) To measure the quasi- 
local density flow every second tunnel coupling was suppressed, 
coupling either odd-even or even-odd pairs, (b) Oscillations of the 
odd-site population in the double-wells with fitted sine-waves for 
t = 100 ^is (solid), 200 /iS (dashed) and 400 /iS (dotted). The value 
of U/J during the relaxation was 5.16(7). (c) Extracted amplitude 
A and phase of the double-well oscillations for odd-even (filled 
circles) and even-odd (open circles) couplings. The solid lines show 
the respective results of the t-DMRG simulations. The dashed lines 
are fits to a linear increase in the phase and a power-law decay of the 
amplitude. The insets show the amplitude in a log-log plot (left) and 
the extracted power-law coefficients (right). The horizontal grey line 
indicates the power-law coefficient a — 0.5 for free and hardcore 
bosons. 



It is key to the experiment that the observed fast damping 
cannot be attributed to a mere classical ensemble averaging 
due to the inhomogeneous distribution of tunnel-couplings in 
the various chains (var( J) / J ~ 0.4%) or the external trap. 
Furthermore, we ensure that the transverse tunnel-coupling 
between adjacent chains J± is always one to two orders of 
magnitude smaller than J. Furthermore, the dynamics of 
a single site - or of the densities of odd sites - cannot be 
described in terms of simple rate equations, and not even 
in terms of Markovian quantum master equations reflecting 
damped motion (see Supplementary Material). Similarly, no 
dynamical mean-field description can capture the dynamics 
for large U ||22| . Hence, any realistic description has to nec- 
essarily include the many-body and non-Markovian features 
of the dynamics, contributing to the challenge for a numerical 
simulation for intermediate times. 

Time-evolution of the quasi-momentum distribution. A 
different view on the relaxation can be obtained from the 
quasi-momentum distribution of the ensemble. When instan- 
taneously switching off all trapping potentials after a relax- 
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FIG. 4. Build-up of short-ranged correlations, (a) Plot of the 
integrated density profiles obtained after ToF versus 'ijt/h for 
U/J = 5.16(7) as obtained in the experiment (left) and recon- 
structed from numerical f-DMRG simulations (right). The images 
show the crossover from a purely Gaussian distribution {t — 0) to 
a more complex quasi-momentum distribution (0 < AJt/h < 2) 
to a purely sinusoidal pattern (AJt/h > 2). (b) Visibility of the 
interference patterns versus ijt/h obtained experimentally (circles) 
and from the simulations (solid curve). The grey line represents the 
measured visibility at AJt/h ~ 5, while the dashed line corresponds 
to the value obtained from the simulation of a homogeneous system 
1 19|. (c) Steady-state value of the visibility measured at AJt/h ~ 5. 
The solid line is a guide for the eye oc J/U . 



ation time t and letting the cloud expand freely for a time 
tToF, the density distribution takes the form nToF(r) c>c 
\wi^{rnY /ht)\'^S{rnY /ht). Here, it;o(k) is the Fourier trans- 
form of the on-site Wannier orbital and the interference term 
for the ensemble of decoupled Hubbard chains in the far- 
field limit is 5(k) = E{jv} E^. y e'^'^i-^'^'\a]ay) with d = 
Xj;s/2 being the lattice spacing along the chain direction. In 
Fig. 4a, we plot the measured density profiles integrated over 
the y- and z-direction as a function of the relaxation time (left 
panel) together with the corresponding patterns reconstructed 
from t-DMRG simulations for the full distribution of chains 
(right panel) for [// J ~ 5. Both the experimental data and the 
numerical calculation show a rapid build-up of short-range co- 
herence, not present in the initial state. 

At short relaxation times AJt/h < 2, the simulation data 
shows a strong cosinusoidal component with a period of 
2hk — 2/itToF / {'m^xs) and weaker contributions from higher 
harmonics. While the former correspond to next-neighbor 
coherences in the system, the latter are a signature of corre- 
spondingly longer range coherences which rapidly decay in 
the relaxation process lfT9l . Due to the noise on the experi- 
mental data, the higher frequency components are weak, but 
can still be identified. For longer relaxation times AJt/h > 2, 



only the next-neighbor coherences remain as also found from 
t-DMRG simulations of homogeneous Hubbard chains with 
finite interactions |[T9l . We extract the visibility of the lowest 
frequency component as described in Ref. 1 14 1 both from the 
experimental data and the i-DMRG calculations (see Fig. 4b). 
In the absence of local currents in the system, the visibility of 
the lowest Fourier component is given by 4E{7v}Re(a|aj_|_j^). 
We find good agreement between experiment and numerics. 
The visibility builds up towards a first maximum at AJt/h ~ 
0.5 corresponding to the first maximum in 7iodd(i) (Fig. 2c), 
followed by dampened oscillations. From AJt/h ~ 2 on, we 
find the visibility to increase slowly with t towards an equi- 
librium value which is significantly higher than the one found 
for a homogeneous system (dashed line in Fig. 4b). In the 
presence of the trap, the system can expand after the quantum 
quench, converting kinetic energy proportional to — (fljaj+i) 
into potential energy associated with {fijp). As a conse- 
quence, the absolute value of (ajflj+i) increases. 

We observe the same qualitative behavior for the whole 
range of interaction strength U / J accessed in the experiments, 
but with a strength of the next-neighbor correlations which 
depends on U / J . In Fig. 4c, we plot the visibility of the in- 
terference patterns found from experiments at different values 
of U/J with at a relaxation time of AJt/h — 5, where den- 
sities and cuiTents are fully relaxed. We find the visibility 
to have a maximum around U / J ~ A while it decreases to- 
wards the analytic limits of U ^ {) and U ^ co. In both 
of these limits, the problem can be mapped onto free particles 
where no coherences will survive for long times. In the lat- 
ter case of hardcore bosons, the coherences are found to be 
suppressed with J/U in the experiments. The same behavior 
can be found from perturbation theory in J/U for a thermal 
state in the lattice |19|. Recent theoretical results show that 
for non-degenerate spectra, the time-averaged state is always 
identical to the maximum entropy state given by the full set of 
constants of motion [23 1. On the other hand, if a state relaxes 
locally, it has to be to the projection of this time-averaged state 
on the respective subsystem. However, it is up to now unclear 
how to physically interpret and identify all constants of mo- 
tion defining the maximum entropy state. Interestingly, the 
findings presented here are compatible with the expectations 
for a Gibbs state defined only by the total energy and the total 
number of particles. Finally note that the build-up of short- 
range phase coherence to a finite value complements the ob- 
served decay of the density pattern and the cuiTents. It cannot 
be explained by any classical dephasing mechanism, but is a 
result of genuine many-body dynamics in the system. 

For the time-evolution of the densities, cuiTents and coher- 
ences, t-DMRG simulations and experimental results show 
a remarkable congruence. This emphasizes the clean im- 
plementation of controlled quantum dynamics in the one- 
dimensional interacting Hubbard-model and the high fidelity 
of the initial-state preparation. All parameters are calculated 
ab inito from the experimental control parameters. Therefore, 
the experiments can be seen as a self-sustained dynamical 
quantum simulation where the simulation effort is the same 
for each value of t and each set of parameters. As long as 
the time-evolution is not perturbed by experimental imper- 
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fection, this dynamical quantum simulator outperforms any 
continuous-time numerical simulation for which the calcula- 
tional effort increases with t. Simulation methods on classical 
computers such as matrix-product state based time-dependent 
DMRG used here suffer from an extensive increase in entan- 
glement entropy which limits the relaxation times accessible 
in the calculations f 24l|25]| . 

Conclusions and Outlook. In conclusion, we have demon- 
strated measurements on the relaxation of a charge-density 
wave of in a strongly correlated one-dimensional Bose gas 
with varying interactions. Using a bichromatic optical su- 
perlattice, we were able to prepare a patterned density state 
with high fidelity and to induce non-equilibrium dynamics 
by rapidly switching on the tunnel-coupling along the chain. 
We could follow the dynamical evolution of the initial state 
in terms of even- and odd-site densities, local currents and 
short-range correlations visible in the quasi-momentum dis- 
tribution. We have compared our measurements to parameter- 
free <-DMRG simulations, finding excellent agreement and 
identifying contributions from next-nearest neighbor hopping. 
All three observables can be seen as local probes of the sys- 
tem and show a rapid relaxation to steady state values as it is 
predicted for Hubbard-type models by a central limit theorem 
1261 l27j . These steady state values are compatible with the 
system globally being in a maximum entropy state. This idea 



can be seen as a dynamical version of Jaynes' principle which 
could recently be substantiated theoretically for local observ- 
ables or reduced states ifTSl |28l |29l and for two-periodic ob- 
servables 1 19 , 22 1 as considered here. After a finite time, the 
closed quantum system cannot be distinguished locally from 
having reached a global Gibbs state under the constraint of a 
set of macroscopic constants of motion set by the initial state 
|15, 26,29,301. 

A direct measurement of global observables to, e.g., iden- 
tify constants of motion of the dynamics is inhibited by the 
ensemble of chains with various particle numbers in the ex- 
perimental realization. This limitation might be overcome by 
either preparing a single chain with fixed length or by select- 
ing a single chain from the ensemble for detection. This will 
open the way to answer what global state is after all reached in 
the evolution, including possible pre-thermalization [31] and 
experimental studies of the eigenstate thermalization hypothe- 
sis f32'-'34l. It is the hope that the present work triggers further 
experimental studies which address theoretically unsolved key 
questions of non-equilibrium dynamics. 
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A/tot (1000) N 

FIG. 5. a Measured atom number on sites with n — 1 versus the total 
atom number. The solid line represents the condition A^tot = Nn^i. 
The arrow indicates the total atom number used in our experiments, 
b Distribution of particle numbers per chain caused by the harmonic 
confinement calculated for our loading sequence. 



SUPPLEMENTARY MATERIAL 
I. LOADING PROCEDURE AND SUDDEN QUENCH 

We started the loading procedure by ramping up the long 
lattice to a depth of 30 i?^' which results in an array of isolated 
two-dimensional quantum gases. The lattice depth is given in 
units of the respective recoil energy ~ / {2rnXf). Sub- 
sequently, we ramped up the two transverse optical lattices 
with wavelengths of Aj,,z — 844 nm to 30 E^-^, crossing the 
Mott insulator transition in each of the two-dimensional gases 
with a maximal filling of one atom per site (see Supplemen- 
tary Material). We then added the short lattice with a depth of 
30 ii'^f * to the long lattice, forming a bichromatic period-two 
superlattice [14). We set the relative phase between the short 
and the long lattice such that every second short-lattice site co- 
incided with a minimum of the long-lattice potential. Finally 
removing the long lattice completely yielded an ensemble of 
one-dimensional density waves | • • • , 1, 0, 1, 0, 1, • • • ) |7, 14| 
with various particle numbers N. 

We veiified the loading of single atoms by bringing each 
atom into a superposition of the states \F = 1, mp = — 1) and 
F — 2,mp = 1), where atom pairs can undergo hyperfine- 
relaxing collisions which would expel them from the lattice. 
We measured the total number of particles with (iVtot) ™d 
without this filtering step (Nn=i), finding that up to a total 
number of 60 x 10^ atoms, no pairs were removed from the 
lattice (see Fig. |5^). For the experiments, we chose a parti- 
cle number of 45 x lO'^ in order to obtain maximally singly 
occupied sites. 



II. DMRG SIMULATION RESPECTING THE GEOMETRY 
OF THE PROBLEM 

After the loading procedure, the occupied sites form an el- 
lipsoid around the center of the trap. Thus, due to the sup- 
pressed transversal hopping, the experimental setup is de- 
scribed by a two-dimensional approximately circular array 
of one-dimensional chains with different particle numbers 
with iVmax = 43. For a realistic theoretical description. 




4Jt/ h 

FIG. 6. DMRG results for the relaxation of the quasi-local densities 
nodd- The solid lines represent the ensemble averaged t-DMRG data 
plotted in Fig. 2a-d. The dashed lines are the simulation results for a 
single chain in the ensemble with A*' = 43 particles. The parameters 
used in the simulations are the same as in Fig. 2. 



DMRG calculations were performed for several particle num- 
bers ranging from 5 to 43 and weighted according to the dis- 
tribution shown in Fig.|5|3. In Fig.|6] we plot the averaged den- 
sities as obtained from DMRG and shown in Fig. 2 together 
with the results for a single chain with the maximal number 
of particles. Both curves differ only slightly on the timescales 
accessible by the simulations. 

In the experiment, the smooth external trap leads to a slow 
expansion of the cloud along the x-direction, rather than a 
sharp reflection of excitations traveling at a velocity cx J in 
the system. Thus, the length of the one-dimensional lattices 
in the DMRG calculations (^ 121) was chosen such that the 
paiticles do not reach the system boundary during the simu- 
lated times. This absence of a sharp reflection also explains, 
together with the averaging over different chains with differ- 
ent particle numbers for which recurrences would happen at 
different times, why no recurrences of the density wave are 
visible in our experiments. 



III. CURRENT MEASUREMENTS 

We measure the quasi-local currents by suppressing every 
second tunnel coupling in 200 /is with the help of the superlat- 
tice, such as to obtain chains of symmetric double-wells (see 
Fig. 3a). For all cuiTent measurements, we chose the super- 
lattice depths to be 40 E!^^ and 4 E!^'', respectively. In this 
configuration, we have a tunnel coupling within the double- 
wells of Jdw — h X 8 kHz and an interaction strength of 
(C//J)dw — 0.2, such that the particles in the double-well 
will practically tunnel independently. In Fig. 3b, we plot the 
relative population on the left side of the double wells n^, as a 
function of the holdtime ^dw and for three different relaxation 
times as they were evaluated for Fig. 3c. Here, we have mea- 
sured the population on the left and right sides of the double 
wells by the same band mapping technique which was used 
for the density measurements. The phase 4> and the amplitude 
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of the oscillations shown in Fig. 3c are extracted from a simple 
sinusoidal fit to the data. 

If we assume the particles in the double-well formed by the 
sites with index j and j + 1 to evolve under the simple single- 
particle Hamiltonian Huw.j — —•^Dwl^jjo-j+i + h-c), we 
find the oscillation amplitude to be 

A, it) = ^j{{h,{t)) +4Im((a,(i)ta^-+iW))', 

(1) 

where t is the relaxation time. From the simultaneous mea- 
surement of the local densities, it is thus possible to recon- 
struct the bare mass current lui{aj{t)^ aj+\ (t)) through the 
barrier between the two sites as a function of the relaxation 
time t. Furthermore, it is evident from Eq. (1), that when- 
ever the mass current vanishes the amplitude is measuring the 
population imbalance between the sites j and j + 1. At these 
points, the quantities {1± A)/2 reproduce the even- and odd- 
site densities, respectively. Therefore they provide two en- 
velopes to 71odd,even(i)- 

For the phase 4>j{t), we find from the same calculation as 
above 



t>j{t) = arctan 



2Im((a,(f)ta^.+i(f)))' 



(n,(t)) 



T'j + l 



it) 



(2) 



We use these two relations to calculate the ensemble-averaged 
amplitude E{jv}^(t) = |E{jv} X]j ^j(^) ^^P(*'/'i(^))l 
phase E{ ^jH^) = arg[E{ jv} J2j ^»(*) exp{i(f)j (t))] from the 
DMRG simulations (solid lines in Fig. 3b). 



IV. DEVIATION FROM MARKOVIAN AND MEAN-FIELD 
DYNAMICS 

In this section, we show that the time evolution cannot be 
described by a Markovian quantum master equation, signify- 
ing the complex relaxation dynamics beyond a situation that 
can be described in terms of rates. That is to say, each con- 
stituent does not see the rest of a chain as a mere bath it is 
weakly coupled to, but intricate memory effects do play a 
role. To see this, it is sufficient to consider the simple case 
of local relaxation dynamics of a single site and the case of 
[7 = 0. We restrict ourselves to the infinite translationally in- 
variant case and an initial state of a density wave as described 
above — any other setting is only more complex and in general 
also non-Markovian. To follow a Markovian time evolution 
of the reduced state of some odd site j, denoted as p, means 
that its time evolution follows a master equation in Lindblad 
form 



dt 



p{t)=t[p{t),h] 



a, 13 



p(<)}h 



where /i is a Hermitian operator that can be different from the 
free Hamiltonian H, and G and are some arbitrary matri- 
ces, reflecting the influence of the rest of the chain. This is 
the most general form of a master equation when the dynam- 
ical map, mapping the initial state |1)(1| — p(0) i-^ p{t). 



reflects Markovian dynamics. Since this is a bosonic free 
model, one can keep track of the evolution by specifying the 
first and second moments alone. The second moment matrix 
of p{t) is given by, with 6 = bj, X = [b + b^)/\/2 and 
P ^i{b^ ~b)/V2, 



lit) 



2{X^){t) {XP)it) + {Px)it) 

{XP){t) + {Px){t) 2{P^){t) 



and the evolution is given by the Gaussian channel ^{t) = 
A{t)^{Q)A^ {t) + B{t). A straightforward calculation gives 

L/2 

with V{t) = e^'*^, as a matrix exponential. 

This evolution then has to be compared with a Markovian 
time evolution, one that satisfies 

A{s+s') = A{s)A{s'), B{s+s') = A{s)B{s')A{s)'^+B{s) 

as would necessarily be true for a dynamical semi-group gen- 
erated by a Gaussian Markovian master equation. A reason- 
able figure of merit of Markovianity of a Gaussian channel, 
is, for a given time t > 0, how different the dynamical map 
is from a Markovian one. There exists more than one reason- 
able way to quantify such a difference, yet all lead to qual- 
itatively identical results. A physically well-defined way of 
quantifying this is a norm difference from the Jamiolkowski 
covariance matrix of the closest Markovian Gaussian channel 
to the given channel at hand ||35]| . A more pragmatic — but 
in the certification of non-Markovianity equally effective — 
means is, for a time f > 0, to simply compute 



inf 



-A{s)A(s')\\ 



(3) 



with the minimization being performed with respect to s, s' > 
such that s + s' = t. This quantity can easily be seen to 
take positive values for intermediate times, giving rise to a 
bound that can be directly related to observable differences on 
states. Hence, no Markovian dynamics — amounting to rate 
equations — can model the dynamics encountered here, and 
memory and genuine quantum many-body effects have to be 
considered. 

The dynamics found here is also not consistent with a 
mean-field picture for values of U significantly different from 
zero. In such a picture, one looks for a time-dependent self- 
consistent solution for 



for L sites in a mean-field Hamiltonian 



Hit) = -J2 bl+ik + blh+i + + i-l)'-'^it))n 

j 

in a variant of the findings of Ref. f2T\. Using Runge-Kutta 
numerical integration one finds that although for short times 
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and small U, the evolution of densities {hj)(t) is quite com- non-triviality of the dynamics, in that a mean field picture can- 
patible with a mean-field picture, one encounters significant not capture the dynamics at hand, and sophisticated i-DMRG 
deviations for larger U and larger times. This again shows the simulations are necessary. 



